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ABSTRACT 

We address the question of an appropriate choice of basis functions for the self- 
consistent field (SCF) method of simulation of the A^-body problem. Our criterion 
is based on a comparison of the orbits found in A'^-body realizations of analytical 
potential-density models of triaxial galaxies, in which the potential is fitted by the 
SCF method using a variety of basis sets, with those of the original models. Our tests 
refer to maximally triaxial Dehnen 7— models for values of 7 in the range ^ 7 ^ 1, 
i.e. from the harmonic core up to the weak cusp limit. When an TV-body realization 
of a model is fitted by the SCF method, the choice of radial basis functions affects 
significantly the way the potential, forces, or derivatives of the forces are reproduced, 
especially in the central regions of the system. We find that this results in serious 
discrepancies in the relative amounts of chaotic versus regular orbits, or in the dis- 
tributions of the Lyapunov characteristic exponents, as found by different basis sets. 
Numerical tests include the Clutton-Brock (1973) and the Hernquist-Ostriker (1992) 
basis sets, as well as a family of numerical basis sets which are 'close' to the Henquist- 
Ostriker basis set (according to a given definition of distance in the space of basis 
functions). The family of numerical basis sets is parametrized in terms of a quantity 
e which appears in the kernel functions of the Sturm-Liouville equation defining each 
basis set. The Hernquist-Ostriker basis set is the £ = member of the family. We 
demonstrate that grid solutions of the Sturm-Liouville equation yielding numerical 
basis sets (Weinberg 1999) introduce large errors in the variational equations of mo- 
tion. We propose a quantum-mechanical method of solution of the Sturm-Liouville 
equation which overcomes these errors. We finally give criteria for a choice of optimal 
value of e and calculate the latter as a function of the value of 7, i.e., of the power-law 
exponent of the radial density profile at the central regions of the galaxy. 

Key words: stellar dynamics ~ methods: iV-body simulations ~ methods: analytical 
- methods: numerical - galaxies: elliptical and lenticular, cD - galaxies: kinematics 
and dynamics 



1 INTRODUCTION 

The 'self-consistent field' (SCF) method of integrating the gravitational A'^-body problem (Clutton-Brock 1972, 1973; Aoki & 
lye 1978; Allen, Palmer & Papaloizou 1990; Hernquist & Ostriker 1992; Earn and Sellwood 1995; Weinberg 1999) has so far 
proved quite useful in the study of particular classes of stellar systems such as isolated galaxies. A powerful characteristic of 
this method is that the gravitational potential, or the force field, at any point of the particles' configuration space are rendered 
by the code in a closed mathematical form, expressed as a series in a suitable set of basis functions. The main task of the code 
is to calculate the coefficients of a series, the resulting potential of which (say V{r, 6, <j)) in spherical coordinates) corresponds, 
via Poisson equation V^<1? = AivGp, to a smooth density field p{r,d,(f)) that is the continuous limit of the mass distribution 
of the N particles. For a known such distribution p(^r,9,(j)), the series coefiicients are given by definite volume integrals of 



* E-mail:ckalapot@phys. uoa.gr (CK); cefthim@academyofathens.gr (CE) 



2 C. Kalapotharakos, C. Efthymiopoulos & N. Voglis 



quantities depending on both p(r, 6, 4>) and on the basis functions. In a iV-body simulation, however, we do not have a priori 
knowledge of the smooth limit p(r, 9, (f>). Thus, we can only obtain Monte Carlo estimates of the values of the coefficients based 
on the particles' positions, i.e., mass elements dm — p{r, 6, 4>)r'^ sin 9drddd<j) are replaced in all the volume integrals by point 
masses (= the particles) and the integrals are approximated by sums over the particles' positions. It follows that the overall 
complexity of the SCF method is 0{N), i.e. linear on the number of particles. Parallelization is straightforward (Hernquist, 
Sigurdsson & Bryan 1995) and applications involving 10*' - 10'^ particles are tractable on present-day computers. 

Owing to its nature, the SCF method is particularly suitable to the purpose of studying the orbital structure or the global 
dynamics of self-consistent models of galaxies. In general, there are three ways of producing smooth estimates of the potential 
of a self-consistent stellar dynamical system: 

a) One chooses an 'ad hoc' analytical model for the potential and then produces a Schwartzchild-type (1979) self-consistent 
model based on a large number of orbits. In this case we actually make no estimate but simply specify the potential in advance. 
Thus the orbital analysis is relatively easy (see Efthymiopoulos, Voglis & Kalapotharakos 2007, section 4 for an indicative list 
of references) since it only requires to integrate orbits in a fixed potential. However, it is well known that the self-consistent 
models constructed by Schwartzchild's method are non-unique, as many different models can be constructed for the same 
system which vary in the velocity distribution. Furthermore, the stability of the Schwartzchild-type models cannot be a priori 
guaranteed (Merritt 1999; Efthymiopoulos et al. 2007). 

b) One makes an iV-body simulation using one numerical method (e.g. direct summation, TREE, or mesh code) up to 
a final time, e.g., when the system reaches equilibrium. Then one obtains a smooth estimate of the potential at the final 
snapshot using a different method, for example, spline interpolation, fitting with polynomial or rational functions, grid or 
the SCF method. When possible, this method is to be preferred over Schwartzchild's method since iV-body equilibria are by 
definition self-consistent and also stable. An important early example was given by Sparke and Sellwood (1987) in the case 
of a fast rotating barred galaxy, while more recent examples, in the case of elliptical galaxies, were given by Contopoulos, 
Efthymiopoulos & Voglis 2000; Jesseit, Naab & Burkert 2005; Muzzio, Carpintero & WachUn 2005 and Muzzio 2006. Attention 
should however be paid on that the use of two different methods to fit the potential, during and after the A'^-body evolution, 
may introduce discrepancies in the resulting analysis of the orbits. The sensitivity of the orbital analysis on the potential 
approximation was studied by Carpintero and Wachlin (2006). These authors found that the relative amount of chaotic versus 
regular orbits in a stationary triaxial model depends significantly on the way chosen to approximate the potential. As shown 
in the sequel (sections 3 and 4), this occurs even if two approximations are similar, e.g., in the case of the SCF method, 
two not very different basis sets. This can be due, for example, to the form of the basis functions, when this does not fit 
the morphology of the galaxy, or to errors introduced in the equations of motion, when a basis set is numerically calculated 
on a grid (as suggested by Weinberg 1999). In any case, such numerical effects affect the orbital analysis in such a way that 
it becomes unclear up to what extent the orbital structure found in a particular system should be attributed to the real 
dynamics of the system or to numerical features of the code that fits the potential. 

c) One integrates the A/'-body system as in (b) and then uses the same scheme (provided by the A''-body code) to fit 
the potential at the final snapshot after the AT-body run (Pfenniger and Frendli 1991, Udry and Martinet 1994, HoUey- 
Bockelmann et al. 2001, 2002; Contopoulos, Voglis & Kalapotharakos 2002; Voglis, Kalapotharakos & Stavropoulos 2002; 
Athanassoula 2002, Kalapotharakos, Voglis & Contopoulos 2004; Shen & Sellwood 2004; Kalapotharakos & Voglis, 2005; 
Voglis, Stavropoulos & Kalapotharakos 2006). This method has internal consistency to a larger degree than (a) or (b). The 
SCF method is an obvious candidate for such simulations. Nevertheless, as already pointed out, the true 'smooth' limit of 
the distribution of the N particles in an actual A-body simulation is not known a priori (the precise meaning of such a limit 
has even been questioned in the literature, Gurzandyan and Savvidy 1986, Kandrup & Sideris 2003). Thus, the question still 
remains on how credible can the reproduction of the orbital structure of an ideally smooth gravitational system be, when 
using an SCF method to fit the potential that is based solely on data of a discretized version of the system, i.e., of the A^-body 
system. A second relevant question is whether there are means that can be devised so that this credibility be maximized. In 
the present paper we deal, precisely, with the above two questions, foetising our investigation, as regards in particular the 
second question, on the appropriate choice of basis functions for the SCF method and on an appropriate method to calculate 
them (when they are not given analytically). 

We propose the following as a numerical test checking the credibility of a particular SCF method (or of any other 
potential fitting method): 1) Choose a pair of smooth potential-density functions p(r), $(r). 2) Create an A^-body realization 
of the distribution p(r) (this can be done by the Monte-Carlo method). 3) Use the SCF code (or any other code) in order to 
calculate a numerical estimate of the smooth potential $jvs(r) from the A^-body system. This we call the response potential. 
4) Finally, compare the orbital structures of the two systems defined by the potentials $(r), $iVB(r). The code can be 
characterized as reliable if the two potentials yield a similar orbital structure. Besides the orbits themselves, a comparison 
may involve various quantities characterizing the orbital dynamics such as, for example, the calculation of the profile of 
the forces, frequency analysis of the orbits, Poincare phase plots, Lyapunov characteristic exponents or other indicators of 
chaos, etc. Our comparisons in the present refer to the central profiles of the forces and to the distributions of the Lyapunov 
characteristic exponents of the orbits. We also give estimates of the impact of Poisson noise, introduced at step (2), on the 
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results. The impact of discreetness was examined in a somewhat different context by Holley-Bockelmann, Weinberg & Katz 



Wc focus on models of elliptical galaxies, in which an important factor affecting the accuracy of the potential fitting is 
the choice of basis functions for the radial part of the potential expansion (Hernquist & Ostriker 1992; Hozumi & Hernquist 
1995). This is related to the fact that the orbits in such systems are sensitive especially on the radial profile of the forces in 
the central parts of the galaxy. On the basis of their central luminosity profiles, the elliptical galaxies are distinguished in two 
groups (Ferarrese et al. 1994; Lauer et al. 1995): a) the 'core' galaxies, with nearly flat or 'shallow' surface brightness profiles, 
corresponding to power- law density profiles p{r) oc r '' with < 7 < 1 (Pridman & Merritt 1997), and b) the 'power- law' 
galaxies in which 7^1, which are further categorized into those with a 'weak cusp' (7 ~ 1) or 'strong cusp' ((7 ~ 2) (Merritt 
& Fridman 1996). In case (a) the force at the centre is equal to zero, while in case (b) the force is finite, for 7 = 1, or infinite, 
for 7 > 1. These differences in the force field at the centre affect directly the regular or chaotic character of the orbits (see 
Efthymiopoulos et al. 2007 for a review), and imply that a SCF method can only be successful if it reproduces correctly the 
central behavior of the forces. In addition, examples are given in which small numerical errors in the potential are amplified in 
the forces, i.e., derivatives of the potential, and even more in the variational equations, i.e., second derivatives of the potential, 
through which the Lyapunov characteristic exponents of the orbits are evaluated. We find that this may lead to an erroneous 
characterization especially of the regular or chaotic character of the orbits. Furthermore, even if an error appears in only a 
small central region, its presence can affect a large number of orbits, in particular box orbits which pass arbitrarily close to the 
centre. Since the box orbits constitute the backbone of many elliptical galaxies, the dynamical implications of differences in 
the central force field are important at least in these galax;ies. This we check by taking as our basic model a maximally triaxial 
Dehnen (1993) 7— model (Merritt & Pridman 1996) for values of 7 in the range ^ 7 ^ 1, i.e., from the limit of a harmonic 
core up to the 'weak cusp' limit. Such cases arc characterized by the presence of many box orbits, and indeed, we find that 
these orbits are strongly affected by the SCF basis set used to obtain the response potential, mainly because of differences 
in the central force field. This contradicts a claim by Hernquist & Ostriker (1992, section 5.2.1) that such differences are 
"probably not significant from a dynamical point of view" . 

In our tests we consider radial basis sets obtained from the literature, i.e., the Hernquist-Ostriker (1992) and the Clutton- 
Brock (1973) basis sets, but also a family of basis sets computed numerically, as proposed by Weinberg (1999). In computing 
the latter, however, we did not use a grid method to solve the associated Sturm-Liouville boundary value problem (Pruess and 
Fulton 1993) because we demonstrate that such methods often result in large errors appearing in the variational equations 
of motion which make use of the second derivatives of a specified basis set (section 3). Instead, we propose a 'quantum- 
mechanical' method of solution of the Sturm-Liouvillo problem which overcomes these errors. This method is applicable when 
the Sturm-Liouville differential equation to be solved is 'close' to another differential equation for which the solution is known 
analytically (the definition of distance of two differential operators is given in section 3) . In our examples below we use the 
'quantum-mechanical' method in order to calculate numerical basis sets which are 'close' to the Hernquist-Ostriker (1992) basis 
set, but improve, however, the representation of the forces at the central parts of the galaxy. The family is parametrized by a 
quantity e which appears in the Sturm-Liouville differential equation (the Hernquist-Ostriker basis set is the e = member 
of the family). We then explore which member of the family better fits, in the N-body realization, the true dynamics of the 
Dehnen model for a particular value of 7. The latter information is given in terms of a function 5(7) specifying, essentially, 
the choice of optimal basis set as a function of the power-law exponent 7 of the central density profile of the galaxy (when 
7^1). This information can be used a priori, i.e., one may choose the optimal basis set for a given iV-body simulation by 
measuring first the value of 7 (from the A^-body data). 

The paper is organized as follows: section 2 presents the general formalism of the SCF method in spherical coordinates, 
following the same notation as in Weinberg (1999), and then refers to the appearance of errors due to several previously 
mentioned sources. Section 3 describes our 'quantum-mechanical' method of determination of numerical basis sets. Section 4 
contains the main results regarding the choice of an optimal basis set following a comparison of the orbits in various Dehnen 
models and in their respective iV-body response models. Section 5 summarizes the main conclusions of the present study. 



2 THE METHOD 

2.1 Basic Formalism of the SCF method 

We start with the basic formalism of the self-consistent field method in the case of spherical coordinates, following the same 
notation as in Weinberg (1999) for cylindrical coordinates. In the SCF approach, a distribution of particles in ordinary space 

is viewed as a Monte Carlo realization of a smooth density field. This field, which is a continuous and differentiable function 
in space and time is given by the integration of the distribution function / with respect to the velocities: 



(2005). 




velocities 




(1) 
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In the sequel we fix the value of the time t and drop this from the arguments of p. In the SCF method the function p(r) is 
expanded in a truncated series of basis functions in coordinates relevant to the shape of the system under study. In the case of 
elliptical galaxies the usual choice are multipole expansions in spherical coordinates. Let pmonopoieir) be the monopole term 
of the multipole expansion of the density and poo{r) a rough estimate that we make of it. Then, we express the monopole 
term as a truncated series in terms of radial basis functions u„oo{r), i.e.: 



poie{r) ^ poo{r) ^ b„oou„oo{r). 



(2) 



The use of a discrete spectrum of basis functions u„oo{r) (labelled by the 'radial quantum number' n) follows from boundary 
conditions imposed to the system (e.g. finite total mass and/or finite size). The coefficients b„oo are unknown and the main 
task of the A'^-body code is to specify their values. Eq.© refiects our expectation that a linear combination of the functions 
Unooir) can fit the residuals of the true monopole term of the real density with respect to our initial estimate poo{r) which acts 
as an envelope in front of the sum in the r.h.s.. In reality, the fitting is efficient if only a small number of terms are needed in 
((2| (the uppermost limit of n, imposed by the A^-body resolution is Umax ~ 0{N^^^) but in practice we use a number of terms 
which is one order of magnitude smaller than this limit (Palmer 1994)). This, on its turn, depends crucially on the initial 
estimate poo{r), which appears not only directly in Eq.([2} but also, as shown below, indirectly, through the Sturm-Liouville 
differential equation which specifies the basis functions u„oo{r)- At any rate, in the same way as for the monopole term, we can 
make in advance some estimate of the profile of the multipole terms of the density by choosing estimate functions pim{r) and 
then fit the residuals of the true multipole terms with respect to the functions pim{r) via series in respective basis functions 
u„im{r). That is, the density is finally written as: 



p{r,6,(f)) = ^ ^ fenim/5;m(T-)u„i™(r)Fi™'(6l,(?i) 



(3) 



l — Q m — — Z n — 



where Yj™ (&,<;/>) are spherical harmonics (e.g. Binney and Tremaine 1987, pp. 655-656). Similarly as for the monopole term, 
the coefficients bnim are unknown while the functions Unim{r) are specified in advance via solutions of a Sturm-Liouville 
differential equation in which the functions pimir) also appear. 

The determination of the basis set of functions Unimir) is done as follows: Repeating for the gravitational potential the 
same procedure as for the density, i.e., selecting in advance some estimate functions <I>i„i(r) for the profiles of the various 
multipole potential terms we also write the potential as 



^{r,e,4>)=Y. £ c„,™$,™(r)u„,™(r)yr(e,' 



(4) 



l — Q m — — l n — 



with unknown coefficients Cnim, and then match Eqs.(|3]) and Q via Poisson equation. This yields finally (in units in which 
G = 1): 



d / 2^2 dUnlm ^ 

dr V dr 



i{i + i)<i>r. 



^ d_ f2 d^lm 

™ dr V dr 



Unlm = — (47rA„Imr $lmPlm)Unlr> 



(5) 



with A„im ~ bnim/cnim- Equatiou ((U, supplemented with appropriate boundary conditions, is a case of the Sturm-Liouville 
eigenvalue problem 



d_ 

dr 



Pim(r-) 



+ qim(r)u = \wi^{r)u 



(6) 



with 



Pim{r) = r <!>;,„ (r) 



qiAr) = l{l + l)<fL(r) - $i™(r)- r 



Whnir) = -4Tvr'^^ijn{r)pim{r). 



d_ / arfJimW 
dr V dr 



(7a) 
(7b) 
(7c) 



From the form of ([5]) it follows that the functions Unim{r) are given by the eigenfunctions of a Sturm-Liouville differential 
operator 



*~'Lm — , 

dr 



d' 



+ qtm{r) 



acting on functions belonging to a Hilbert space with the inner product definition 

< f\g >= / f{r)g{r)wim{r)dr 



(8) 



(9) 
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where ra,ri, are two radii at which boundary conditions must be given in the form (Pruess & Fulton 1993) 

am — a2 ( pi,n{r)^] = X„i,n ( a'lU — a'i^] at r = Va (10a) 



dr J \ ar 

biu + b2 ^pim{r)^^ ^ at r = rt (10b) 

with constants ai, a2, a'i,a'2,bi,b2. The problem admits a discrete set of solutions, i.e., a discrete set of eigenvalues Xnim and 
eigenvectors Unim, n — 0,1, 2, ... of Cim- In this way, the basis functions Unim{r) are specified by the initial choice of estimate 
functions pimir), ^im{r), which are hereafter called the kernel functions of the Sturm-Liouville problem, and by the boundary 
conditions. The index n is called the radial quantum number. 

In galactic problems, the radii ra,rb are usually set equal to ra = (centre of the system), and rb = Rp or rb oo, 
depending on whether we consider a system ending at a finite radius Rp or at infinity. In the former case, the boundary 
conditions at rb — Rp represent the request of continuity, and of continuous derivative, of the potential function at the point 
Rp where we pass from Poisson to Laplace equation. Such boundary conditions can always be cast in the form (|10|l . 

If the kernel functions $im{r), pim{r) satisfy Poisson equation \7^^im{r) = AirGpimir), then they are called a potential- 
density pair of functions. In that case we always have Aozm — 1 and uohn = const. Independently of whether the kernel functions 
are potential-density pairs or not, the eigenfunctions Unim.{r) of the Sturm-Liouville problem ([5]) are always orthogonal with 
respect to the inner product definition ((9]). 

In an A'^-body simulation we use the above formalism in order to obtain a smooth potential function as follows: 

a) We specify in advance the sets pim{r), $tm{r), and then Unim{r),\nim, as explained above. 

b) Given the particles' positions, we find estimates of the coefficients b„im of an underlying 'smooth' density field, by 
exploiting the orthogonality relation < Unim\un>im >= 3„,„', together with the orthogonality of the spherical harmonic 
functions. Namely, equation ((3} is inverted, yielding the value of any particular coefficient bnim as 

bnim^ / / R{r,e,4>)drded4> (11) 

Jra Jo Jo 

where 

i?(r,e,0) = -~4n'S>i^{r)u„,^{r)Yne,(l>)p{r,9,(l>y sind. 

Assuming now that the positions of the N particles provide a discrete realization of the smooth density field p{r,6,(p), the 
volume integral (|lip can be evaluated by the Monte Carlo method, i.e., as a sum over the particles' positions: 

JV 

6„i™ ~^-'i-K^i^{ri)u„i^{ri)Yr''(Oi,<t)i). (12) 
1=1 

c) We finally calculate the coefficients of the potential series as Cnim ~ bnim/^nim- 

The use of a discrete sum instead of a continuous integration implies that, contrary to what the term 'smooth field' might 
suggest, there is always some numerical noise in the system introduced by discreteness effects, which produces 'relaxation' 
effects in the A'^-body simulation (see Palmer 1994, or Weinberg 1996 for a detailed discussion). In our numerical examples in 
the sequel we examine in detail the effect of this type of noise on the orbits. At any rate, as already mentioned the goodness 
of the fit of the smooth potential of a system by the SCF series depends crucially on the choice of basis functions Unim which, 
on their turn, depend on the choice of appropriate kernel functions pimir), ^imir). To this we now turn our attention. 



2.2 The choice of basis functions 

In order to determine a suitable basis set, a common strategy in the literature is to use kernel functions pimir), <^im{r) such 
that the solutions of ([5]) are reduced to some simple form given, e.g., in terms of special functions or in a closed polynomial 
form. We mention the following examples (see Hernquist & Ostriker 1992 for a detailed review): 
(i) the Clutton-Brock (1973, hereafter CB) set. The kernel functions are: 



pUr) = ^{2l + l)[2l + 3) (13a) 



with boundary conditions 



^ and fim ^^^^^ = 0. 

„ r->oo dr 

r = 



dr 

In the form (|10|) these conditions are given e.g. as ai = a'l = a'2 — 0, a2 — 1, and 61 = 0, 62 = 1. The zeroth order term of the 
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density expansion 

3 a 
" ^(a2 + r2)5/2 

yields a half-mass radius equal to unity, ri/2 = 1, if a — 0.77 ~ 10/13. The corresponding potential is 

1 



The above set pooO)'3?ooo is the well known Plummer model. The eigenfunctions u„jm(r) are Gegenbauer polynomials of the 

^2 _ ^2 

form Ci^^(0 (independent of m) where ^ — 



+ 

(ii) the set of Allen et al. (1990). In this case there are no kernel functions of the potential or density, i.e. 

pim{r) = l, $i„(r) = -l. (14) 

The inner boundary condition is dUnim{Q)/dr = while the outer boundary condition corresponds to the matching of the 
solutions of the Poisson and Laplace equations at some finite radius = Ftp 



dUnlmir) 



dr 



= 0, r +{l+ l)Mnim(r) 

ar 



= 0. 

r=Rr, 



r=0 

In the form (|10p we set ai = a'l = a'2 = 0, a2 — 1, and bi = I + 1, b2 — 1/Rp. The eigenfunctions Unimir) are spherical Bessel 
functions. This set better fits systems with a harmonic core. 

(iii) the Hernquist-Ostriker (1992, hereafter HO) set. The kernel functions are 

, , /—I (2; + l)(; + l) ar' .^^ . 

Pim{r) = V47r — '-^ '- 15a 

2tt r (o + r)^'+'^ 



(a + r) 

The zeroth order term of the density 



'^Ur) = -V4^ ,„ , ^,21+1 - (15b) 



Pooo{r) = \/47r-!--!- 



27r r (a + r)^ 

yields a half-mass radius equal to unity, ri/2 = 1, if a = (1 \/2)'^^'^ which is our choice of value of a in the sequel. The 

boundary conditions are = finite at = and the same as in the CB set at infinity. In the form (llOp we give 

dr 

2Z + — 

precisely the same constants as in the CB set. The eigenfunctions Unim are Gegenbauer polynomials of the form Cn ^ (0 

where f — -. 

r + a 

The choice of an appropriate set of basis functions is related to the morphological characteristics of the system under 



study. For example, if the kernel function poo(f') is selected as in the HO basis set (Ea. (ll5a|l 'l. this function yields a power-law 
cusp poo{r) ~ r~^at the centre, so the HO set cannot easily fit systems in which the density cusp at the centre is shallower 
than this law. In fact, Hernquist & Ostriker (1992, subsection 2.3) give an idealized example in which a flat profile at the 
centre can be represented by a linear combination of the n = and n = 1 monopole terms of the HO basis set, because the 
coefficients of these terms are balanced in a way so as to eliminate the power-law dependence r~^ . However, we can show that 
if, due, for example, to the Poisson noise, the balance is slightly distorted in an actual A^-body calculation, the singular 
behavior of the density at the centre reappears as a result of the numerical error. In our normalization units and constants 
(a = (1 + \/2)~^^'^, which corresponds to a = 1 in Hernquist & Ostriker 1992), the example refers to the spherical density 
profile: 

Po(?-) = ^ (16) 
47r (a + r)* 

which is clearly flat at the centre. While the kernel function poo{r) given by Ea. (|15a[l goes as poo{r) ~ r~^ at the centre, a 
combination of the two first monopole terms of the radial basis set allows one to remove the singularity and obtain the profile 
(|16p. Precisely, we have: 



_ a'+i2^'+^r(2;+|) / nl{n + 2l+l) ^21+i (r- 

""""" 4^ Y (l + 20(/ + l)r(n-f4/ + 3)^" \r + a^ 

so that 

P''' = 2^r(a + r)3 

and 

^a3/2(r - a) 
2iTr(a + r)* 
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Thus 



PoooW + cpioo(r) 



27rr(a + r)'^ 



c\/5(r — a) 



and, ii c — l/%/5 (which corresponds to c = 1/12 in the units of Hernquist & Ostriker (1992)), we find 

1 



pooo + ^Pioo = ^(^^ oc p„(r). 



(18) 



However, the balance in Ea. (|18|l is quite sensitive to numerical errors. For example, if, due to discreteness effects, there is a 
small error eg in the ratio c of the numerical coefficients of pooo and pioo, i-e., c = 1/y/E + eo, the error introduced in the 
evaluation of the density is 



1 , VSa^/^ , 1 I 3£o(r-a) 

pooo + (—1= + eoJPlOO = PQ(r) = —, ; + - —. ' 



(19) 



and we have limr-,o Po(^) = i.e. the density as calculated by the SCF method with the HO basis set becomes singular at 
the centre. 

Fig. 1 shows a numerical calculation of the above effect. Using a (50 x 100 x 200) spherical polar grid in (6^, (j), r) respectively, 
logarithmic in the radii (from r = to r = 15), and linear in the angles, we produced an A'^-body realization of the spherical 
system (|16p with A'' — 1.25 x 10^ particles as follows: a) a theoretical value of the number of particles in each grid cell 
was calculated, that is given by ANth ~ po{rc)r^sm0cArA6A(f), where rc,6c are the values of r and 6 at the centre of the 
cell, b) A random integer number AA^ is then determined which has a Poisson distribution corresponding to a mean value 
< AN >= ANth and dispersion a = AN^l'^ . c) The values AA'^ for all the cells were normalized so that the total number 
of particles is kept equal to A'^. d) Each cell was filled with AA'^ particles located randomly within the cell with a uniform 
distribution. Finally, we used this distribution of particles to calculate the coefficients of the ten first monopole terms &nOO, 
n = 0, . . . , 9, of the sum ((2]) with the HO basis set. 

The first two coefficients in the above simulation obtain values yielding a ratio 6000/&100 = 2.27381, which is close to, but 
not exactly equal to the theoretical ratio feooo/^'ioo ~ "s/S — 2.236067. In this case the relative error of the two coefficients is 
due to the fact that the A'^-body system is truncated at a radius r = 15 after which the density field (|16p . for the used number 
of particles and grid, is seriously undersampled (e.g. we have less than 1 particles per cell). This error in the sampling can 
be partly reduced by introducing variable masses of the particles (Sigurdsson, Hernquist & Quinlan 1995) so that there are 
more particles tracing the smooth mass distribution at the outer radii. Even so, however, the error cannot be reduced further 
from the Poisson noise limit caused by the discrete number of particles. This is estimated as follows: We evaluate the integral 



fEa. (|ll|l by truncating the radius at rj, — 15 (instead of infinity) and find theoretical coefficients foooo^rjiTs') 



^theoretical 

Oioo,r6=i5 irom 



-2.5x10' 



-5.x 10' 



-7.5 X 10' 



-l.xIO 




Figure 1. (a) The profile of the radial force Fr(r) corresponding to a mass distribution with the density profile I II6I I (with 1.25 X 10*^ 
particles) as derived analytically (solid curve) and via the SCF calculation by the first two monopole functions (n = 0, 1, I = m = 0) 
(dashed red curve) and by the first ten monopole functions {n = — 9, I = m = 0) (dotted blue curve) of the HO radial basis set 
implemented on an A^-body realization of this mass distribution (see text for details), (b) A zoom of (a) in the central region of the 
model. For large radii the force is represented very good in both cases (two and ten monopole terms expansion). For small radii the first 
two monopole terms yield a small but definitely finite value of the force at r = while the first ten monopole terms yield a conspicuous 
error in the force near the centre r = 0. 
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the theoretical density p{r, 8, 



po{r), and also numerical coefficients by the sums (I12|) yielding a relative error 



Ab, 



'000,rb = 15 



^theoretical 
"000,ri, = 15 



2.7 X 10" 



Ab 



100,rj=15 



^theoretical 
"l00,rb = 15 



8.4 X 10" 



Both these relative errors are comparable to a Poisson noise error of order 0{1/VN), with iV ~ 10^ 

Now, both errors due to the truncation and to the Poisson noise contribute in that, while the overall fitting of the radial 
force profile by the HO basis set with the two first terms (n = 0, 1) is quite satisfactory (Fig. la), the force obtained at the 
centre by the same terms is finite (Fig. lb). Even so, one may argue that the finite value of the force is rather small. However, 
this only happens because we considered the first two terms in the radial expansion, while, in a typical A''-body simulation 
in which we are not aware of whether the density profile is close to some idealized example, we typically use many more 
terms (say up to n = 9). Theoretically, the coefficients of these terms in Eq. (|18p are equal to zero. In practice, however, we 
find that these terms also have small non-zero values which, because of the 1/r singularity in the HO kernel function poo{r), 
increase dramatically the error in the central force (Fig. lb). Furthermore, while the appearance of the error can in principle 
be reduced to very small values of r, it always afi^ects orbits which pass arbitrarily close to the centre, such as the box orbits. 
This, and other numerical effects, will be examined in the subsequent sections. 



3 A 'QUANTUM MECHANICAL' METHOD OF DETERMINATION OF NUMERICAL BASIS SETS 

Weinberg (1999) stressed the need for an adaptive algorithm producing basis sets tailored to the morphological details of 
the system under study, so that the series representation of the potential has good convergence properties and numerical 
instabilities such as that mentioned in the previous section are avoided. Translated to spherical coordinates, Weinberg's 
proposal is essentially the following: 

a) Choose kernel functions pim{r),^im{r) which are as good estimators as possible of the profiles of the A^'-body system 
to be run. 

b) Solve numerically the Sturm-Liouville boundary value problem (|6]) by special solvers such as the SLEDGE code (Pruess 
& Fulton 1993). Such solvers are based on variants of the shooting method, and they yield the solutions u„im in tabulated 
form, i.e., at the points of a grid. 

In principle, Weinberg's proposal extends considerably the applicability of the SCF method by enlarging the freedom of 
choice of kernel functions ^imir), pim{r), which can even be adaptive, i.e., change in the course of an A^'-body simulation. 
However, we show now that the use of tabulated values on a grid imposes restrictions to accuracy so that we have to devise 
an alternative method for the solution of the Sturm-Liouville problem. 

To this end, we point out first that when the functions Unim are known only on a grid, the first and second derivatives of 
the potential can be calculated directly only by a finite difference scheme. However, the use of any such scheme whatsoever 
cancels the property of smoothness of the potential. In order to retain the latter we can use an interpolating function for 
each tabulated function u„im, with sufficient number of continuous derivatives, that joins the values of any function Unim{r) 
at successive points of the radial grid. However, if there are small numerical errors at the grid points, these errors also show 
up when one calculates the derivatives either by finite differencing or by the derivatives of the interpolating function. Fig. 2 
shows how does the error appear when using interpolation. The SLEDGE solver was used in order to produce numerically 
the HO basis set Unim, which is given analytically by Ea. (|17p with I = 0, ... ,4 and n = 0, . . . , 9, so that comparisons between 
the numerical and analytical solution can be made. The solution was required on a grid of 2.2 x 10'* points in the interval 
^ r ^ 22, with a tolerance 10~*. After the basis functions were calculated, each function was interpolated by a cubic spline, 
so as to secure the continuity of the first and second derivative of the interpolating function at the grid points. Fig. 2a shows 
the interpolated numerical solution in the case n = 1 and in the interval 0.1 ^ r ^ 0.2. Although SLEDGE could not reach 
the requested input tolerance 10~** in many parts of the solution, the error actually rendered by SLEDGE in the values of 
the function mioo was still relatively small (of order 10"'' or below. Fig. 2d,e). However, this small error is amplified when one 

calculates the derivatives (Fig. 2b) via the interpolating function. Then the error becomes, in general, of order 10~^, 

reaching up to 10~^ at some particular points of the numerical solution (Fig. 2d,e). In fact, a careful examination of these 
points showed that the grid solution had small jumps there, of order 10~^, and these were forcing the interpolating function 
to introduce artificial points at which the convexity was seriously distorted (a function could even turn locally from convex 

to concave). Thus, when calculating the second derivative 5 — , the error was dramatically amplified at these particular 

points (Fig. 2c). We see that the error locally becomes of order unity, while the overall error in the second derivatives ranges 
from 10~* to 10~^, i.e., similarly as for the forces. We found similar numerical effects in practically all the functions u„i„i that 
were determined numerically by the SLEDGE solver. Furthermore, when orbits were calculated in an TV-body realization of 
a maximally triaxial Dehnen 7— model, for 7 — 0.4 (see next section), using the HO basis set calculated numerically up to 
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Figure 2. (a) The HO eigenfunction Mioo(f') (n = 1,1 = m = 0) versus r as derived by tlie SLEDGE solver of the corresponding 
Sturm-Liouville problem, after a cubic spline interpolation of the numerical values on a grid of stepsize dr = 0.001. (b) The derivative of 
the curve (a) as calculated by the interpolating function. Two non smooth spots are just visible at r ~ 0.11 and r ~ 0.175 respectively, (c) 
Second derivative of (a). The error at the two spots is now clearly amplified, creating large errors in the variational equations of motion 
in which the second derivatives of the eigenfunction appear, (d), (e) The absolute and relative errors of the interpolated numerical 
solution with respect to the analytical solution in logarithmic scale for the function «ioo (black solid line), and for the derivatives 
duioQ{r)/dr (red dotted line) and d'^u{r)ioQ(r) / dr'^ (blue dashed line). The error in the first derivative is amplified with respect to the 
error in the eigenfunction, while the error in the second derivative becomes of order unity at the spots, (f ) The evolution of the finite 
time Lyapunov characteristic number (LCN) for one orbit in the maximally triaxial Dehnen model with 7 = 0.4 (see section 4), initial 
conditions xq ~ 0.045760, 3/0 — 0.004794, ~ 0.001598, VxO = 0, Vyo = 0, v^o = 0, when the potential is reproduced by the HO basis 
set calculated analytically (blue solid curve) or numerically (by the SLEDGE code) on a grid (red dashed curve). The expansions in both 
cases reach up to Umax = 9 and Imax = 4. The LCN curve shows the orbit to be regular in the first case, and chaotic in the second case 
with LCN » 10-2-5, 



timax ~ 9, Imax = 4, we found orblts that were affected by these errors to such an extent that artificial positive Lyapunov 
characteristic numbers were introduced in cases of orbits that were in fact regular (Fig. 2f). 

In order to overcome the above numerical difflculties, we suggest now a 'quantum-mechanical' method of solution of 
the Sturm-Liouville problem ([6} which is similar to the quantum-mechanical perturbation theory of bound eigenstates (e.g. 
Merzbacher, 1961, pp. 413-437). This method is applicable when a numerical basis set is required that is not very different 
from another basis set known analytically . Let u^i^{r) be the functions of the known set. These are solutions of the problem 
([6| for a Sturm-Liouville operator jCjJ^' which is determined by a choice of kernel functions as, e.g., in subsection 2.2, and 
for specific boundary conditions of the form (|10[l . Let, now, dm be a different Sturm-Liouville operator corresponding to a 
different choice of kernel functions. We seek to solve the eigenfunction problem 

£-imUim{r) = \u'imir)uim{r) (20) 

i.e., find the spectrum of successive eigenvalues X„i and eigenvectors Un'im{f) of 'C-im, n' = 0, 1, for which we request to 
satisfy the same boundary conditions as those satisfied by the functions u^^i^{r). Since the latter form a complete basis of 
the space of functions with the given boundary conditions, any function can be written as a linear combination of them. We 
thus write 

oo 

U„'lmir) = J2 d-r^',nlmU^°Li'r)- (21) 
n = 
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In view of ((20]), Eq.((2T]) takes the form 

OO OO 

dn',nim/;;m«i°L('^) = ^n' X] ,r^lmWln,{r)u^°irn{r) ■ (22) 
n=0 n=0 

Multiplying both sides of (I22II with — — j-^ , integrating over all radii from ra to r%, and recalling the orthogonality of 

the functions u'^i^ir) yields the linear set of equations: 

00 

X] < ^fL\-Mlrn\u'^L > d„',nlm = K' d„' Jim (23) 
n=() 

where the linear operator A4im is defined as 

Mim = ! jCim (24) 

Wim(r-) 

and 

,(0) 1 A/(. i„(o) ^= P „(") r.^A/(, „(o) 



< -Z1>1<"^I<"L / (25) 
The problem (|25|) is equivalent to the original eigenvalue problem and it can be written in the matrix form: 

Ti.im ■ d„/^;„ = Xd„i.lm (26) 
where d„/ is a column vector with entries equal to the coefficients d„tni,^, that is d^/ — [dri' fiim,dn' ,iim, d„t 2i„i, ■ ■ ■), and 
Him is a matrix with entries Ti.im,ij =< ''J-um\-^im\ufim ^- This formally corresponds to what is referred to as the 'Hamiltonian 
matrix' in quantum mechanics, while Ea. (|26|l is analogous to the quantum mechanical procedure of diagonalization of the 
Hamiltonian matrix. 

The numerical steps in order to solve the Sturm-Liouville problem (I20|l are then summarized as follows: 

a) Starting from a known basis set u^^i^) calculate the integrals (|25p and hence the matrix Him- In fact, the matrix Him 
contains an infinite number of entries, thus we can only compute k x k truncations of Him- Nevertheless, the determination 
of the eigenvalues and eigenvectors of (I26|l for different truncation orders is a convergent procedure (numerical examples 
are given below). The convergence is faster when the set of eigenfunctions u^^^^ and u^im correspond to 'nearby' operators 
Ml°2,Mim. By this we mean that the kernel functions (r), "I>im(''), and w,'^ (r), «;;m(r), through which the operators 
■^hn. a'^'^ Mim are defined, should have a small distance in their functional space, the latter being defined for two arbitrary 
functions /, g as, for example, the euclidian distance 

rf'(/,5)= r{f{r)-g{r)fdr. (27) 

b) Solve the eigenvalue problem (|26|) for the kxk truncated problem. This determines eigenvalues X^nm and eigenvectors 
d„',im. The latter are translated to the new eigenfunctions u^iimir) = d^,;^ ■ UjJ^(r), where U;J^(r) is the column matrix 




Figure 3. (a) The numerical value of the eigenvalue A900 of the new basis set as a function of the rank k (order of truncation of the 
matrix Ti) when e = 0.02. Note that we need at least a 10 X 10 matrix in order to start calculating this eigenvalue, (b) A zoom of (a) 
in which the convergence of the eigenvalue A900 to a final value, as k increases, is shown in greater detail. The final value is practically 
reached when fc > 20. 
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with the functions u^^^^{r),n — 0, 1,2, ...fc — 1 as entries. Notice that the matrix A4im is real and symmetric, thus its 
diagonahzation is a numericaUy fast and accurate procedure. 

We have implemented the above procedure in order to produce a numerical basis set that satisfies the following two 
properties: 

i) It is 'nearby' to the HO basis set in the sense of small functional distance of the kernel functions given by Eq. (|27p . 

ii) It can reproduce density profiles which are shallower in the centre than p(r) oc . 

The new basis set was defined as follows: Selecting the initial kernel functions "I>;^ (r) , (r) and basis functions u^^i^ 
to be the HO set, we define the new kernel functions 

■f;™M = <°iW (28a) 
r— 1 (21 + l)(l + 1) ar' , , , 

Thus, the only change with respect to the HO set is the introduction of a softening parameter e in the singular factor l/r of 
the HO density kernel function. It follows that the operators Cl'2^,Cim are identical ZI;^ = Cim, but the operators mI'2_, Mim 
are different, Mf^ ^ Mim, because 'wl'^{r) / wim{r). Furthermore, 

lim d{wl2,wi,„) = 

thus, for s sufficiently small, the operators A1;^, Alim are 'nearby' according to the previously given definition. 

Fig. 3 shows the convergence of the numerical solution of the eigenvalue problem (|26|l for the above kernel functions, as 
a function of the truncation order k. The abscissa in Fig. 3a,b is the dimension k of the k x k matrix TiYm^ truncation rpj^^ 
ordinate shows the numerical value Agoo (in the case e — 0.02) of the eigenvalue of the tenth eigenvector of Ti^m^ truncation 
which corresponds to the quantum numbers n — 9, and I — m — 0. Clearly, when k is as small as fc = 11, the numerical value 
of Agoo calculated from the matrix Tigo^'^ truncation jg already very close to the value Agoo — 62.29235 corresponding to the 
limit k — > oo. A zoom of Fig. 3a is shown in Fig. 3b, in which we see that the limiting value of Aggg is practically reached 
after k > 20. In fact, we always find that the limiting values of the eigenvalues from Xoim up to Xnim are specified essentially 
up to the computer's double precision limit when one uses matrices TCim truncated at an order k — 3n. 

Fig. 4 shows an example of numerical basis function {u52m{r) for n — 5, I = 2, independent of m) calculated for different 
values of the parameter e, compared to the HO and CB basis functions for the same quantum numbers. The plot focuses on 
the region of inner radii (r ^ 0.25). Clearly, when r > 0.15 all the numerical basis functions are much closer to the HO basis 
set than to the CB basis set, while, at radii below r = 0.15 the effect of introducing e is manifested, namely a basis function 
becomes less steep at the centre as the value of e increases. Thus, by choosing different values of e we can better control 
different behaviors of the central potential or density profiles of a simulated system. We also note that the CB basis function 
deviates considerably at the centre from both the HO and the numerical basis functions, a fact expected since the CB kernel 
function (Plummer sphere) represents systems which are rather flat at the centre. 




Figure 4. The forms of the radial basis functions 1152™, (f) (n = 5, i = 2, independent of m) in the inner region (r < 0.25) in the HO 
and CB cases and in the case of the 'quantum-mechanical' calculation for various values of the parameter e. 
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4 A NUMERICAL TEST. THE DEGREE OF ORDER AND CHAOS VIA THE SCF METHOD 

In the present section our task is to test the accuracy of reproduction of the orbital content of a model elliptical galaxy by 
various SCF codes differing in the choice of radial basis set. The model is Denhen's (1993) 7— model with an ellipsoidal radius 
m (Merritt 1999, section 1). The density reads: 

p (^) = ^^~^\^°^ m-^(r„ + m)-(*-^) (29) 
inabc 

with ^ 7 < 3, where 

= 4 + + 4 (30) 
a-^ 

is the ellipsoidal radius corresponding to a triaxial system with axial ratios a : b : c, a ^ b ^ c > Q, and M equal to the 
total mass of the galaxy. The parameter 7 determines the exponent of the power-law profile of the density at the centre, i.e., 
p{m) ~ m ' at the centre which essentially yields also a radial profile p(r) ~ r~^. We are interested in the case of 'core' 
galaxies in which 7 is in the range 0^7^ 1. The potential corresponding to the density (|29p reads (Merritt & Fridman 
1996): 



2(2 - 7)r„ ^ Jo y'{T + a2)(r + b^)(T + c^) 

with 



<^{x,y,z)^~ ^^^ X ^ ,,..L,.....,.\. dr (31) 



In order to produce an A'^-body realization of the previous system, we work as in the numerical example of subsection 
(2.2) using = 1.25 x 10^ particles arranged in a (50 x 100 x 200), spherical polar grid. We found that such a number of 
particles was necessary because the results regarding all the tests below were becoming robust against the number of particles 
for A'^ of the order of A*' = 10^ or higher. This fact is related to various effects caused by the Poisson noise (subsection 4.2 
below). The unit of length is determined so that the half-mass radius is at ri/2 = 1, while the A''-body system is truncated 
at a radius r = 15, containing 95% of the total mass of the model system (in which the distribution of the mass extends 
theoretically to the limit r 00). 

In all the models we set a = 1, b = 0.7905, c = 0.5 corresponding to a maximally triaxial model. Six different Dehnen 
models of progressively higher power-law exponents were examined, namely 7 = (perfectly harmonic core), 7 = 0.2, 0.4, 0.6, 
0.8, and 1 (weak cusp). After the A'^-body realization for each of these models was produced, the potential and density were 
fitted by the SCF method using eight different radial basis sets. These are the HO and the CB sets, as well as the basis sets 
derived by the 'quantum-mechanical' method of section 3, for the values of the parameter e equal to e = 0.005, 0.01, 0.02, 
0.04, 0.06 and 0.08. The latter are called e— modified basis sets (modified with respect to the HO basis set). 
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Figure 5. (a) The x axis component Fx of the force along the direction x = y = z = r for the Dehnen model 7 = 0.4 (black solid 
line) and the approximations of the same force component obtained by implementing the SCF method to the Af-body realization of the 
Dehnen model using the HO, CB and the e-modified basis sets (the values of e are indicated in the figure. The number of particles in 
the Af-body realization are 1.25 X 10'^. (b) a focus of (a) near the centre. The force with HO has a finite value at r = 0, while the force 
with CB is zero at the centre but its overall central profile is much smoother than in reality. On the other hand, all the e-modified basis 
sets behave better than the HO or CB sets at the centre. The optimal value of e by visual inspection is close to e = 0.02. 



Appropriate SCF basis sets for orbital studies of galaxies 13 



There are two different numerical tests performed in these systems: a) we compare the accuracy of reproduction of the 
behavior of the forces of the Dehnen model at the centre, for each choice of SCF basis set, and b) we compute a library 
of 1200 orbits and compare the number of orbits that are found to be regular or chaotic, as well as the distribution of 
the Lyapunov characteristic numbers produced by the numerical integration of the variational equations of motion in each 
potential representation. In the sequel we separately analyze these two categories of numerical tests and compare their outcome 
as regards the 'optimal' basis set to be used, i.e., the value of e for which the Dehnen model and SCF agreement are better. 
This information is given as a function of the value of 7. 
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Figure 6. The quantity ^ versus e for each 7— model examined. The quantity AF^ is the mean square difference, inside 

a — (x.y,z) 

a central sphere of radius r = 0.01, between the force component a (where a can be x, y, z) of the specific 7— model and the same force 
component as derived from different SCF basis sets implemented to the TV-body realization of the Dehnen model. The quantity F^-i^ is 
a characteristic value of the force for all the models set equal to Fch = 10^. The red and the blue horizontal lines correspond to the HO 
and CB basis sets, respectively, (a) 7 = 0, (b) 7 = 0.2, (c) 7 = 0.4, (d) 7 = 0.6, (e) 7 = 0.8, (f) 7 = 1. For 7 = the minimum error 
is for £ = 0.04. As 7 increases the optimal value of e decreases. After 7 = 0.8, the HO basis set gives smaller error than the e— modified 
basis sets, which, however, remains in relatively higher levels than the minimum error found for smaller values of 7. The CB basis set 
gives the worst results (as regards the force error) for all the values of 7. 
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4.1 Forces 

Fig. 5 shows the reproduction of the forces by the eight different SCF basis sets in the A'^-body reahzation of the 7 = 0.4 model, 
compared to the analytically derived forces through the differentiation of (|3ip with respect to the coordinate variables. In all 
calculations below the ten first radial basis functions n = 0, 1, ... 9 are used for each angular function , with Z = 0, . . . 4, 
m — —I, . . . ,1. We found that despite the large increase in the number of basis functions there was practically no significant 
difference observed when n was raised up to Umax = 13 and I up to Imax ~ 6. The odd angular functions Z = 1, 3 or m = 1, 3 
are kept in the SCF calculation despite the fact that theoretically the coefficients of the corresponding terms in the potential 
or density expansion are zero (because the Dehnen model is fully symmetric with respect to the three principal planes (x — y), 
{x ~ z), OT {y — z)). Numerically, these coefficients turn to have small values that yield an estimate of the effect of the Poisson 
noise of the A''-body realization on the numerical values of the coefficients. In Fig. 5, The F^, projection of the force is plotted 
as a function of r when the force is calculated in the direction x — y = z = r. Fig. 5a refers to a radial extent up to r = 1, i.e., 
up to the half-mass radius of the system. We see that all the SCF models provide in general a satisfactory reproduction of the 
true force in this range, except for the CB model for which the agreement is good only in the outer radii (r > 0.2). On the 
other hand, the various models are diversified in their behavior close to the centre, which is seen by zooming to ^ r ^5 0.1 in 
Fig. 5b. We notice immediately the failure of both the HO and CB sets to provide a reasonable representation of the forces in 
this scale. The force, as derived by the HO set, turns to have a significantly non-zero value at the centre. This is precisely the 
problem of the non-perfect balancing of the coefficients analyzed in subsection (2.2). On the other hand, the CB set yields a 
zero force at the centre, but near the centre the force is seriously softened with respect to the true force. One obtains a better 
agreement if one uses the e— modified basis sets, as exemplified by the plots of Fig. 5b for e = 0.01, 0.02, or 0.06. 

In order to quantify the quality of the fit of the forces by the different SCF models, we cannot use as a relevant 
quantity the fractional errors AF{r)/F{r), as a function of r, because we have F{r) as r ^ 0, yielding fractional errors 
AF{r) /F{r) ^ 00 in this limit. We thus use as a relevant index the quantity 

where AFJ is the average value of the squared differences of the forces in the axis a = x,y, z in a, spherical volume ^ r ^ 0.01 
and Fch is a characteristic value of the force set equal to Fch = 10^, used to normalize the errors of all the models. The 
quantity < -^J- > for the eight SCF models, and for all the different 7— models is shown in Fig. 6 (the case 7 = 0.4 of 
Fig. 5 corresponds to the panel 6c). In the ideal 'harmonic core' case (Fig. 6a, 7 = 0), a significantly non-zero value of e 
{^optimal — 0.04) is required in order to minimize the error in the central forces to a level < >~ 10"'^. On the other 

hand, as 7 increases, the optimal value SopUmai, at which the force error is minimized, is shifted towards smaller values of e 
{Eoptimai — 0.02 when 7 = 0.2 (Fig. 6b), Eopumai — 0.01 when 7 = 0.4 (Fig. 6c), EopUmai — 0.005 when 7 = 0.6 (Fig. 6d)). 
When 7=1, the e— modified basis set yields worse results than the HO basis set, while it is still better than the CB set 
(Fig. 6e,f). A typical estimate of force errors with the HO basis set is < -^3- >~ 10~^. We emphasize that such errors in the 
force determination appear inside a very small region of the galaxy (up to r = 0.01 when the half- mass radius is r = 1). Thus, 
one may claim that the errors are not important in the overall A^-body simulation. Nevertheless, these errors are important 
when one calculates the regular or chaotic character of the orbits, as demonstrated in the next subsection. 



4.2 Lyapunov exponents and the regular or chaotic character of the orbits 

In order to check numerically the regular or chaotic character of the orbits, we create, for each 7— model, a library of 1200 
orbits calculated by a set of initial conditions uniformly distributed on four different equipotential surfaces of the 7— model 
potential ^^{x,y,z) corresponding to the values of the energy Ei = 0.95"I>7 (0, 0, 0), E2 = O.S"!?^ (0, 0, 0), E3 = 0.6<1?t, (0, 0, 0) 
and E4 — 0.4$^ (0, 0, 0). In the energies Ei and E2, which are close to the central value of the potential, we find many 'box' 
orbits which are regular. The existence of box orbits is actually guaranteed by the fact that the force at the centre is zero for 
7 < 1. On the other hand, as the value of the energy increases the phase space is dominated by different families of tube orbits 
or chaotic orbits. The tube orbits follow essentially the classification of Statler (1987) for the perfect ellipsoid. In general, 
chaos increases as the value of 7 increases, a fact mainly associated with the destruction of the regular character of the box 
orbits. 

In order to characterize the orbits as regular or chaotic we use the same numerical criterion as in Voglis et al. (2002). 
This is a combination of two numerical indices, called the 'finite-time-Lyapunov number' Lj and the 'Alignment index' Alj 
respectively (see Voglis et al. 2002 for details). The numerical values of these indices are calculated for each orbit, labelled by 
the index j, over an integration time equal to T = 1200 radial periods. In a two-dimensional plot of log(yl/j) versus log(Lj) 
(Fig. 7, for 7 = 0.6), the regular orbits accumulate in the right and middle part of the diagram, in an almost straight horizontal 
segment around the value log(Lj) ~ log(l/r) ~ —3. On the contrary, the chaotic orbits occupy the up and left part of the 
diagram and have a considerable scatter in the values of Lj, which are a measure of the maximal Lyapunov characteristic 
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Figure 7. The method used to distinguish between regular and chaotic orbits. The Alignment Index Alj, j = 1, 1200 of each orbit 
is plotted against the finite time Lyapunov number Lj in log — log scale. A detailed description of this method is given in Voglis et al. 
(2002). The specific plot refers to the case of Dehnen model with 7 = 0.6. Each point in this plot corresponds to an orbit that was 
integrated for 1200 radial periods. Regular orbits have low values of Lj (decreasing with the integration time as and high values of 
Alj, thus the group of regular lies in a region near log(Lj) ~ —3 and log(j4Jj) > —3.5. Chaotic orbits have large values of Lj and very 
small values of Alj (Alj < 10"^"). The orbits on a lane joining the above two groups are weakly chaotic orbits (i.e. sticky orbits) which 
have just started exhibiting their asymptotic chaotic behavior. 



exponent of the orbits. In between these two groups of points there is a lane of points connecting the two groups. This refers 
to weakly chaotic, or 'sticky' orbits, and as T increases, there is a continuous, albeit decreasing, flow of points through the 

lane towards the group of chaotic orbits. However, this flow is irrelevant in the characterization of the orbits for times greater 
than T, because it should imply that some orbits characterized as 'chaotic' have, in fact, Lyapunov times much longer than 
the age of the galaxy. 

Fig. 8 shows, now, the main result as regards the choice of an appropriate set of basis functions that better reproduces the 
regular or chaotic character of the orbits. All the panels show the distributions of the logarithm of the flnite-time Lyapunov 
numbers Lj in the case 7 = 0.6 for the orbits which are within or above the transport lane of Fig. 7, i.e., the orbits characterized 
as chaotic. The distributions are not normalized, i.e., the area below each distribution is equal to the percentage of orbits that 
were characterized as chaotic. Thus, differences in the total area covered by two distributions mean a difference in the total 
percentage of the orbits characterized as chaotic, while differences in the shape reflect a different normalized distribution of 
Lyapunov exponents. The distribution shown with solid line refers to the precise calculation in the analytic Dehnen force field, 
while dashed plots refer to computations with different SCF basis sets. Furthermore, in this figure all the SCF distributions 
refer to a calculation in which we switched off the odd terms of the angular expansion of the potential as derived by the SCF 
code for reasons explained immediately below. 

The following are some basic remarks regarding these diagrams: 

a) The computation with the HO basis set overestimates by a factor of about 1.7 the percentage of chaotic orbits in this 
model. Thus, as shown below, while the true percentage for 7 = 0.6 is ~ 22%, the HO fit yields ~ 38%. In fact, the difference 
in the two percentages is even higher if one takes into account the odd terms of the angular expansion of the potential (in this 
case we find a percentage 48% with the HO basis set). However, we know that the true value of these terms in the Dehnen 
model is equal to zero, so that non-zero values can only be associated with asymmetries in the number of particles on the 
two sides of any of the principal planes of the galaxy induced by the Poisson noise of the A''-body realization. The simplest 
way to measure the latter is by measuring the size of the coefficients of the odd angular functions / = l,3orm = l,3in the 
potential expansion. We thus use the following measure of the Poisson noise: 

_ ||coefficients of odd angular terms|| 
1 1 all coefficients! I 

where the norm || • || means sum of the absolute values of the coefficients. In the case of Fig. 8a, we find PN ~ 2.3 x 10~^ 

which is consistent with an estimate PN ^ 0{l/\fN) with A'^ = 1.25 x lO'^. This looks like a relatively small number, which 
however produced a 10% difference in the percentage of chaos with the HO simulation. 

b) The computation with the CB set (Fig. 8h) underestimates the percentage of chaotic orbits by a factor ~ 2.75 (22% 
true percentage against 8% with CB). The underestimate is even more serious (6%) if one switches off the odd angular terms. 
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Figure 8. The distribution of the values of log(L) of the chaotic orbits for the Dehnen model 7 = 0.6 (same black solid line in each 
panel) compared to those found by the HO, the e— modified (for the indicated values of e) and the CB basis sets (blue dashed lines) in 
the A'^-body realization of the same system. From these figures it is clear that HO overestimates chaos, CB underestimates chaos and the 
best estimation is obtained by an e— modified basis set for £ between e = 0.005 and £ = 0.01. Odd angular terms m = 1, 3 are excluded 
from the computation because they are affected by Poisson- noise asymmetries in the A'^-body realization (see text for details). 
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Figure 9. The percentages of chaotic orbits in the original Dehnen models and with the various SCF basis sets applied to the cor- 
responding A'^-body realizations. Each panel indicates the corresponding value of 7. Black solid lines correspond to the percentages of 
chaotic orbits in the original Dehnen models, pale blue (dotted) and blue (dashed) lines to the calculation with CB, with and without 
odd angular terms, orange (dotted-dashed) and red (dashed) lines to the calculation with HO, with and without the odd angular terms. 
Black solid curves correspond to the calculation with the £— modified basis sets for the various values of e (curves with open circles = 
calculation with the full angular expansion, curves with solid rectangles = calculation without odd angular terms). The optimal value 
of e with respect to the criterion of agreement in the percentage of chaos is determined by the point where the curve of the e— modified 
calculations (curve with solid rectangles) intersects the horizontal black line (original Dehnen model) 
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Figure 10. The percentages of chaotic orbits found by the use of the optimal basis set (see Fig. 9), for each value of 7, when the 
non-symmetric (odd) angular terms in the potential expansion axe turned on (lower curve) or off (upper curve). 

The forms of the two distributions are also quite different, a fact meaning that the CB fit is rather unsuitable to represent a 
galaxy with central power-law exponent 7 = 0.6 (or beyond). 

c) The best results are found when we use an s— modified basis set with eopUmai between the values SopUmai = 0.005 and 
Soptimai = 0.01 (Fig. 8b,c). For this value of e both the percentage of chaotic orbits and the distribution of the Lyapunov 
exponents derived by the SCF method are in good agreement with the true percentage and distribution. 

Fig. 9 shows a comparison of the true percentages of chaotic orbits with those found by the SCF fit with different basis 
sets, for all the examined values of 7. An optimal value of e, based on the 'percentage of chaos' criterion, is determined 
by the point where the curves with open dots intersect the horizontal solid line, marking the true percentage. This value, 
however, is also contaminated by chaos due to the Poisson noise on the odd angular terms, and a better determination is 
made when these terms are switched off (solid curves with black rectangles). In any case. Fig. 9 renders immediately clear 
that the HO fit produces overestimates of the percentage of chaotic orbits in the whole range of values ^ 7 ^ 1, while the 
CB fit underestimates the percentage of chaos for the values 0.4 ^ 7. Up to 7 = 0.8 an optimal £— modified model can be 
found yielding the best agreement with the true percentage. Near this value of 7, however, the situation is reversed, and in 
the limit 7 ^ 1 it is the e— modified models yielding underestimates and the HO model yielding the best results. This is, 
precisely what is expected by noticing the power- law profile of the HO monopole terms at the centre. 

Fig. 10 shows the comparison of the percentages of chaotic orbits found by the use of the optimal basis set, for each 
value of 7, when the non-symmetric (odd) angular terms in the potential expansion are turned on or off. In all cases, the 
numerical noise due to the non-symmetric (odd) terms increases the percentage of chaotic orbits. This fact justifies the use 
of symmetrized N-Body realizations of a system, as e.g. by HoUey-Bockelmann et al. (2001), when an orbital analysis is 
requested. 

Fig. 11 shows a comparison of the distributions of the finite-time Lyapunov exponents of the chaotic orbits for all the 

7— models considered, with the distributions derived via the HO model, optimal e— modified, and CB model. This figure checks 
essentially whether, at the value of e at which the agreement of the percentages of chaotic orbits is good, the agreement of 
the distributions of the finite time Lyapunov numbers is also good. This check is necessary, because it is possible that two 
very different distributions yet cover the same total area. The plots in the middle column of Fig. 11 clearly show that the 
distributions with the £— modified basis set, at the optimal value e = SopUmai, are indeed in good agreement with the true 
distributions, except in the case 7 = 1 in which the best agreement is obtained by the HO basis set. 

In order to check the robustness of the above results on the particular way chosen to create the N-Body realization of 
a Dehnen 7— model, we created a different realization, in which there is no use of polar grid, but an acceptance-rejection 
Monte Carlo algorithm was utilized to produce an N-Body system for the 7 = 0.4 model. The mass contained in a small 
volume element corresponding to values of the spatial coordinates in the intervals (m, m -I- dm), {6, 6 + d9), and (0, (f> + dcp) is 
given by dm = m^ahcp{m) sm{9)dm dO d(j), implying that the mass distribution is separable and uniform in the coordinates 
Xi = m^/3, X2 = — cos(^) and X3 = 0. We thus readily obtain acceptance-rejection Monte Carlo realizations of the models. 
In practice we transform m to ^ = (m — l)/(m -|- 1) so that when ^ is given values in —1 ^ < 1, m can take arbitrarily 
large values. In reality, given the finite number of particles in the simulation, the distribution of the matter is undersampled 
at large radii. We thus checked the robustness of our results against three different truncation radii, namely a) rt = 15 (which 
is the standard choice in the previous simulations), b) n = 30, and c) rt = 100. The results are summarized in Fig. 12a, 
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Figure 11. The log{L) distributions of the chaotic orbits as derived by the HO, optimal e— modified, and CB approximations compared 
to the true distribution of the corresponding Dehnen model. Each row refers to one indicated value of 7. Black solid curves correspond 
to the real distribution while blue dashed curves are the distributions with HO (left column), optimal £— modified (middle column), or 
CB (right column) we have for the specific value of 7. For 7 ^ 0.6 the agreement is better between the real and optimal e— modified 
calculations. For 7 = 0.8 the results with e— modified {e = 0.005) and HO (e = 0) are comparable, while for 7 = 1.0 the best result is 
with HO. 
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Figure 12. (a) The log(L) distributions of the chaotic orbits, as derived by the HO basis set (with all terms), in a rejection-acceptance 
Monte Carlo realization of the 7 = 0.4 Dehnen model, when the truncation radius is rt = 15 (red solid line), rt = 30 (blue dotted 
line), rt = 100 (green dashed line). The real distribution by the original 7-model is plotted by black solid line, (b) The absolute error 
\I{rt) — I{oo)\ versus rt in log — log scale. The red solid line corresponds to evaluations of I{rt) using Monte Carlo realizations with 
TV = 10^ particles while the blue dashed line to Monte Carlo realizations with N = lO"^ particles. We see that beyond a value of rt the 
error due to Poisson noise dominates the error due to finite radius rt- 




Figure 13. (a) The optimal value of e versus 7 based on the 'forces' criterion (see Fig. 6). (b) the optimal value of e versus 7 based on 
the criterion of the percentage and log(L) distribution of the chaotic orbits (see Figs. 9 and 10). 



for the distributions of the finite-time Lyapunov numbers calculated with the HO basis set. When the truncation radius is 
at rt = 30, there is a marginal improvement of the results (the total percentage of chaotic orbits found is 37%, against 41% 
when rt = 15, and 14% real percentage). But even with a much higher truncation radius {rt = 100) the improvement is still 
small (32% against 14% real percentage of chaotic orbits). We thus conclude that the distributions shown in Fig. 11 do not 
change appreciably by changing the sampling technique or the truncation radius. 

A simple argument can show why the scaling of the numerical error with the truncation radius rt saturates at some 
value of rt. This is based on equating the error in the evaluation of generalized integrals, due to a truncation of the limits 
of integration, with the error due to the Poisson noise at large radii by a Monte Carlo evaluation of the integral. Namely, 
since the Monte Carlo sampling uses a finite number of particles, the drop of the density at large radii implies that the mass 
distribution is seriously undersampled at these radii. An example is given in Fig. 12b, which refers to a calculation of the 
integral 

lin) = / '47rr'p(r)$oo(r)dr, (34) 
Jo 

corresponding essentially to a truncation at finite radius of the generalized integral (|ll|l for a spherical Dehnen-model p(r) — 
r '(1 -I- r)~''* when the oscillatory behavior of the functions u„oo{r) is neglected. The integral (|34|l was evaluated by a 
Monte Carlo sum using TV = 10^ or A'^ = 10^ particles. Fig. 12b shows the absolute error \I{rt) — ^(00) 1 when a different 
Monte Carlo realization of the system is produced for each value of rt- For small rt, the error is large and it is due to the 
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truncation of the generalized integral at insufficiently small radius. However, beyond a value of rt, the Poisson noise due to 
the undersampling of the density at the outer parts clearly dominates, causing variations of the error which are of the same 
order as the error due to the finite truncation. When A'^ = 10^, the critical value of rt is estimated as rt = 8 (in units of the 
half-mass radius), and the error beyond this radius fluctuates between 10~^ and lO^"*. On the other hand, for iV = 10^, the 
truncation radius beyond which the Poisson noise dominates raises to rt = 25, and the error beyond this radius fluctuates 
between 10~^ and 10~^. Those values arc compatible with the estimates given independently in section 2. 

Finally, Fig. 13 shows a comparison of the optimal value of e determined in the previous analysis by (a) the central force 
criterion (subsection 3.1), and (b) the chaotic percentage criterion (present subsection). The latter criterion suggests the use 
of somewhat larger values of e than by the former criterion, when 7 < 0.4, while the two criteria yield nearly equal values 
for 7 > 0.4. At any rate, the main conclusion drawn from the above numerical examples is that the characterization of the 
regular or chaotic character of the orbits of an A'^-body system depends crucially on the accurate numerical representation 
of the 'smooth' potential that presumably underlies the mass distribution of the N particles. In the framework of the SCF 
method, this conclusion is translated into the need for very careful choice of basis set. In the above examples we used a 
'quantum-mechanical' method in order to calculate numerical basis sets which are close to the HO basis set, but they can 
better fit the behavior of all smooth quantities (potential, forces, density) at the centre when the power-law central density 
profile is shallower than r^^. Such profiles naturally arise in Af-body simulations of the remnants of galaxy mergers (e.g. 
Jesseit et al. 2005 and references there in). Nevertheless, the method is applicable, with different starting basis sets, to a much 
wider class of stellar dynamical systems and the cautions raised in the present section are, very probably, relevant to such 
systems as well. 



5 CONCLUSIONS 

In the present paper we address the question of the choice of an optimal basis set of functions for orbital studies of galaxies 
simulated via the self-consistent field (SCF) method. Our approach is to create A/'-body realizations of a Dehnen analytical 
7— model of a triaxial galaxy and then reproduce the gravitational potential by the SCF method, using different basis sets 
which are either analytical, i.e., the Hernquist-Ostriker and Clutton-Brock sets, or numerical, depending on a parameter e 
that modifies the behavior of the potential at the centre with respect to the HO basis set. We then compare a number of 
quantities characterizing the orbits (central force profiles, percentage of chaotic versus regular orbits and distributions of the 
Lyapunov characteristic numbers) in the original model and in the SCF reproduction of the potential. Our main conclusions 
are the following: 

1) If the kernel functions (section 2) used in the construction of a radial basis set have a singular behavior at the centre, 
this basis set is unsuitable for simulating galaxies which do not exhibit the same singular behavior. In particular, even if the 
singularity can be theoretically removed by a balancing of the coefficients of some terms in the SCF series, the numerical 
noise induced by the Monte-Carlo evaluation of the values of the coefficients destroys the balance and results in large errors 
appearing in the evaluation of the forces or derivatives of the forces, especially in the central parts of the galaxy. 

2) When following the methodology suggested by Weinberg (1999) for the determination of numerical radial basis sets, 
shooting methods yielding tabulated values of the basis functions on a grid should be avoided, because they, too, introduce 
largo errors in the evaluation of the derivatives of the forces, yielding large inaccuracies in the integration of the variational 
equations of motion. 

3) We propose a 'quantum-mechanical' method of determination of numerical radial basis sets that overcomes the above 

problems. This method is based on the diagonalization of a truncated matrix which results from a spectral formulation of 
the Sturm-Liouville problem, equivalent to the procedure referred to in quantum mechanics as the diagonalization of the 
Hamiltonian matrix. The stability and accuracy of this scheme is demonstrated by speciflc numerical examples. 

4) We study the orbits in a family of triaxial Dohnon 7— models for various values of 7 in the range 0^7^ 1, corresponding 
to the case of elliptical galaxies with 'shallow' central cusps. For each 7— model we use a family of different radial basis functions 
which are modifications of the HO basis set, derived by implementing the 'quantum- mechanical' algorithm of (3). The basis 
functions depend on one small parameter e defined so as to yield the HO basis set in the limit e ^ 0. The numerical tests 
concern comparisons of a) the behavior of forces in the central parts of the galaxy, b) the percentage of regular and chaotic 
orbits, and c) the distribution of the Lyapunov exponents of the chaotic orbits, between the exact 7— model and the various 
A'^-body - SCF realizations using the HO, CB and e— modified models for various values of e. When the standard basis sets 
(HO or CB) are used, we find large deviations in all three criteria between the exact and SCF results, except for the HO basis 
set for 7 = 1. On the other hand, the best agreement is found by using e— modified basis sets for values of e in the range 
0.005 ^ e ^ 0.06 (in units of the half mass radius). The optimal value of e is given as a function of the value of 7. The latter 
can in principle be determined in advance (before or during the A^-body simulation) by measuring the power-law exponent of 
the central radial profile of the system. 

All the numerical codes used in the present paper are available by the authors upon request. 
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